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, The mass spectrum of heavy pseudoscalar mesons, described as quark-antiquark 

cn 

' bound systems, is considered within the Bethe-Salpeter formahsm with momentum- 



dependent masses of the constituents. This dependence is found by solving the 
Schwinger-Dyson equation for quark propagators in rainbow-ladder approximation. 
Such an approximation is known to provide both a fast convergence of numerical 

X" 

^ ' methods and accurate results for lightest mesons. However, as the meson mass 



increases, the method becomes less stable and special attention must be devoted to 
details of numerical means of solving the corresponding equations. We focus on the 
pseudoscalar sector and show that our numerical scheme describes fairly accurately 
the vr, K, D, Dg and rjc ground states. Excited states are considered as well. Our 
calculations are directly related to the future physics programme at FAIR. 
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I. INTRODUCTION 

The description of mesons as quark-antiquark bound states within the framework of 
the Bethe-Salpeter equation with momentum dependent quark masses, determined bv the 
Schwinger-Dyson equation, is able to explain successfully many spectroscopic data [l|-l9|. 
Contrarily to traditional phenomenological models, like quark bag models, the presented for- 
malism maintains important features of QCD, such as dynamical chiral symmetry breaking, 
dynamical quark dressing, requirements of the renormalization group theory etc. The main 
ingredients here are the full quark-gluon vertex function and the dressed gluon propagator, 
the calculation of which is entirely determined by the running coupling and the bare quark 
masses. In principle, if one were able to solve the Schwinger-Dyson equation within all pQCD 
orders, the approach would not depend on any free parameters. However, due to known 
technical problems, one restricts oneself to calculations of the few first terms of the perturba- 
tive series, usually up to the one-loop approximation. The obtained results, which formally 
obey all the fundamental requirements of the theory, are then considered as exact ones with, 
however, effective parameters. This is known as the rainbow-ladder approximation for the 
Schwinger-Dyson equation. The merit of the approach is that, once the effective parameters 
are fixed, the whole spectrum of known mesons is supposed to be described on the same 
footing. Also the excited states are supposed to be described within this approximation. 

It should be noted that there exists other approaches based on the same physical ideas 
but not so sophisticated, e.g. employing simpler interactions, such as a separable interaction 
for the effective coupling p]. Such approaches describe fairly well the properties of light 
mesons, nevertheless, investigation of heavier mesons and excited states, consisting even of 
light (u,d,s) quarks, requires implementations of more accurate numerical methods to solve 



the corresponding equations. Among other successful efforts in this realm the Refs. 10|, 111 ] 
must be also mentioned. 

In the present note we are going to apply the combined Schwinger-Dyson and Bethe- 
Salpeter (BS) formalisms to describe the meson mass spectrum including heavy mesons and 
excited states as well. Particular attention is paid to the charm sector which, together with 
the baryon spectrosco py, is a major topic in the FAIR research programme. Two large collab- 



orations at FAIR 
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13| plan precision measurements. Note, that it becomes now possible 
to experimentally investigate not only the mass spectrum of the mentioned mesons, but also 
different processes of their decay, which are directly connected with fundamental QCD prob- 
lems (e.g., U{1) axial anomaly, transition form factors etc.) and with the known problem of 
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changing the meson characteristics at finite temperatures. The latter is crucial in understand- 
ing the di-lepton yields in nucleus-nucleus collisions at, e.g. HADES. All these circumstances 
require an adequate theoretical foundation to describe the meson spectrum and the meson 
covariant wave functions (i.e. the BS partial amplitudes) needed in calculations of the relevant 
Feynman diagrams and observables. 



II. PROPAGATORS AND SCHWINGER-DYSON EQUATION 

The Bethe-Salpeter and Schwinger-Dyson equations in Minkowski space contain poles and 
branch-point singularities which strongly hinder the procedure of finding numerical solutions. 
Usually, to avoid these difficulties, one performs the Wick rotation and formulates the corre- 
sponding equations in Euclidean space, where all singularities in amplitudes and propagators 
are removed, so that the equations can be solved numerically. The known Mandelstam tech- 
nique allows then to calculate matrix elements of observables which, being analytical functions 
of the relative energy, are the same in both Minkowski and Euclidean spaces. 

In our case we consider the Schwinger-Dyson equation for the quark propagator within 
pQCD with summing all diagrams up to one-loop. In calculations of diagrams the chiral 



symmetry breaking is implemented from the very beginning joj. The exact results, even 
only up to one-loop diagrams, after proper regularization and renormalization procedures are 
rather cumbersome for further numerical calculations. Nevertheless, in practical calculations 
one can employ reasonable parametrizations for the corresponding vertices and propagators 
to find solutions numerically. So, a phenomenological expression, inspired by calculations 
of the mentioned diagrams and preserving the requirements of the theory, for the combined 
running coupling and gluon propagator has been suggested in Ref. ^ 



ln[r + \ 



^QCD 



where the first term originates from the infrared (IR) part of the interaction determined by 
non-pertubative effects, while the second one ensures the correct ultraviolet asymptotics. Ac- 
cordingly to the known fact that the contribution of the IR part is predominant for formation 
of bound states, in what follows we neglect the second term and restrict ourselves to the IR 
one. Then, as seen from Eq. ([T]), we are left with only two effective parameters, D and u. 

The calculation of the renormalized Feynman diagrams leads to a fermion propagator de- 
pending on two additional parameters. In canonical calculations these are the renormalization 
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s-'(p) = n-pA(p) + B(py, s,(p) = ^^l™^^ . (2) 



constant Z2 and the self-energy corrections Usually, for further simplifications of calcu- 

lations, instead of Z2 and one introduces other two quantities A{p) and B{p) in terms 
of which the quark propagator Sq{p) reads as 

-27 ■ pA{p) + B{p) 
p'^A^ip) + B^{p) 

Then with such a representation of the quark propagator the Schwinger- Dyson equation 
in Euclidean space has the form (cf. Refs. |7|, l9|) 

S;\p)=z^.p + m + ^ I ^^[g'Dip-l)]^^^,S,{lh,, (3) 

where rfi is the bare quark mass and the effective kernel D{p ~ I) is 

D(k\„ = ^^e-"''"' (S,. - . (4) 

Note that ([3]) is a four dimensional integral equation. To solve it one usually decomposes the 
kernel over a complete set of basis functions, performs analyticaly some angular integrations 
and considers a new system of equations relative to such a partial decomposition. In our 
calculations we expand the interaction kernel into hyperspherical harmonics 

D{{p-iy) = D(J9^/^H|/|cos7) = 5^D„(j9,/)C^(cos7) , (5) 

n 

where are the Chebyshev polynomials and D'^{p,l) denotes the corresponding part of 

the kernel D{k'^)/k'^. Eventually, the system of equations to be solved reads as 
1 4 /• dU^ A{1) ^ , 

/ + - ^^^(^'0 + 5Z.;(p,/))] , 

f dip Bil) ^ , „ 

The resulting system of equations on ([6]) is a system of one-dimensional integrals and can 
be solved numerically, e.g. by an iteration method. We found that iteration procedure for ([6]) 
converges rather fast and practically does not depend up on the choice of the trial functions 
for A{p) and -B(p). 



III. BETHE-SALPETER VERTEX FUNCTION 

To determine the bound state mass of a quark- ant iquark pair one needs to solve the Bethe- 
Salpeter equation. In rainbow-ladder approximation and in Euclidean space it reads [?, 9| 

r(P,p) = {-\) I ^,i,S{k^)mk)S{k^h^{g'D{p-k)),^, (7) 
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with r being the Bethe-Salpeter vertex function. The color structure has been factorized 
exphcitly. P means the total momentum of the bound state, p{k) is the relative momentum 
within the quark pair and k+ = k + C,P, k_ = k + — 1)P. The result does not depend on ^ 
due to the covariance of the Bethe-Salpeter equation. 

Equation is written in matrix form, i.e. the vertex function T{P, k) is a 4x4 matrix and, 
therefore, may contain 16 different functions. The general structure of the vertex functions 
describing bound states of spinor particles has been investigated in detail, for example, in 
[l^ . r can be expanded into spin - angular momentum functions: 

r(P4,p) = ^ga{p4Ap\)Ta{p) (8) 

a 

and in functions ga{p4, \p\) in its turn the dependence on the angle x can be written as 

gM, \p\) = ^gl{p)Xj,o{i){xp)- (9) 

j 

This allows to perform the angular integration leaving us with a set of coupled one- dimensional 
linear integral equations for the functions g^{p) 

9fip) = l i^E4'"(P'%r(^) (10) 

where p{k) is the length of the Euclidean 4-momentum, p = \/pi^ + p^. The matrices A^p{p, k) 
can be calculated analytically. 

The series ([9]) converges rather fast and in practice only few terms need to be taken into 
account. Then by choosing a suitable method of integration, e.g. Gaussian quadrature, 
the system of equations ( ITOl) can be written as a system of homogeneous linear equations. 
Schematically, it can be written in the form 

g = Kg. (11) 

The condition det(i^' — 1) = is sufficient for the existence of a bound state. Hence, formally 
the zeros of the determinant (lei{K — 1) determine the solution of the BS equation, including 
also excited states. 

Here an important moment is worth to be emphasized. Usually, when solving the BS 
equation for constituent particles, i.e. for particles with constant masses, the resulting system 
of partial equations is real. In case of momentum dependent masses the Bethe-Salpeter 
equation becomes complex and requires the knowledge of the quark propagator for complex 
momenta k± (where k± = ^P±k are the momenta of quarks) and, hence, requires to solve the 
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Schwinger-Dyson equation also for complex momenta. For small meson masses (e.g. M ^ 500 
MeV) only a small energy region contributes to the integral in Eq. f lTOj) . k < 1 GeV. In this 
case the propagator functions can be obtained by using the solution for real momenta, and 
afterwards A{p) and B{p) are calculated at complex momenta from ([6]) For heavier 

states, M > 1 GeV, the imaginary part of the quark momenta, Im k±, becomes rather 
large and the integrand in ([6]) rapidly oscillates as a function of fc, hindering an accurate 
computation of the integral. 

To avoid this problem a witty trick has been suggested in Ref. Q]. It is based on the 
observation that, as seen from the Schwinger-Dyson equation, all the values of the momenta 
/c^ are located within a domain limited by a parabola. 

Then, due to the Cauchy's theorem, it suffices to know the solution^ along this parabola to 
be able to compute it everywhere inside the corresponding domain 



IV. NUMERICAL RESULTS 

As an example of our numerical study we exhibit in Fig. [1] the energy of the lowest bound 
states of a hypothetical meson qq^ consisting of one given quark q with the mass known from 
the Schwinger-Dyson equation, bound with a second quark q^ for which the input bare mass 
fhx is let to vary arbitrarily. The corresponding effective parameters have been chosen as 
uj = 0.5 GeV and D = \Q GeV~ and the bare masses for q correspond to (m, rf, s) 

quarks, q=u (with = 0.005 GeV), q=s (with rhs = 0.115 GeV) and q=c (with rfic = 1.0 
GeV). This figure illustrates the whole mass spectrum of pseudoscalar mesons with masses 
up to 3 GeV/c^. So, if the q^ quark corresponds to a c quark, then at the intersection of 
the vertical line = rric 1 GeV/c^) with "ux" curve one obtains the D -meson (with 
the quark contents auc), with the " sx" curve the Dg meson and with the "cx" curve - the rjc 
meson, respectively. It is worth noting that the "ux" curve crosses the m„ = 0.115 GeV line 
roughly at the same value of Mgg^ as the " sx" curve crosses the rhu = 0.005 GeV line, thus 
proving a check of consistency of the approach and, at the same time, describing correctly 
the lowest pseudoscalar us state corresponding to the K meson. It can be seen that even 
without a fine tuning of ^ the meson mass spectrum is reproduced fairly well: 135 MeV 
(7r° meson), 497 MeV {K meson), 1870 MeV (D^ meson), 1970 MeV {Df meson) and 2980 
MeV {rjc meson). 

Clearly, a special parameterization of the gluon propagators with two adjusted quantities 
and the self consistent determination of three bare quark masses serve as input for obtaining 
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FIG. 1: The bound state masses of a system qq^ as a function of the bare mass rhx, the other quark 
being either a u quark with rhu = 0.005 GeV (the curve labeled as "ux") or a strange quark with 
rhs = 0.115 GeV (the curve labeled as "sx") or a charm quark with rhc = 1.0 GeV (the curve labeled 
as "cx"). 

finally five pseudoscalar meson masses. The straight forward extension to the scalar, vector 
and axialvector mesons without further adjustments will provide the ultimate test for the sub- 
tleties of the numerical implementation of the employed theoretical scheme and the inherent 
approximations. Work along this line is in progress. 

V. CONCLUSION 

The method of solving the Schwinger-Dyson equation in rainbow-ladder approximation 
in Euclidean space by using the hyperspherical harmonics basis is proposed. The obtained 
numerical solutions are then used to solve the Bethe-Salpeter equation for the meson mass 
spectrum in a large interval of meson masses. In solving the Bethe-Salpeter equation a new 
set of basis functions has been used which allows a further easy decomposition of the Bethe- 
Salpeter vertex functions into hyperspherical harmonics basis 

The obtained mass spectrum for pseudoscalar mesons in a wide range, ranging from pions 
to rjc mesons, is in a good agreement with experimental data. Excited states were considered as 
well and found to be also in a good agreement with experimental data and with calculations 
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by other groups. By solving the Bethe-Salpeter equation, the corresponding partial wave 
functions are also obtained which will allow, in future, to calculate a variety of observables 
related to physical programmes at, e.g. FAIR. 
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